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ABSTRACT 


The  description  of  a  beam  element  by  only  the  displacement  of  its  center  line  leads  to  some 
difficulties  in  the  representation  of  the  torsion  and  shear  effects.  For  instance  such  a  rep¬ 
resentation  does  not  capture  the  rotation  of  the  beam  as  a  rigid  body  about  its  own  axis. 
This  problem  was  circumvented  in  the  literature  by  using  a  local  coordinate  system  in  the 
incremental  finite  element  method  or  by  using  the  multibody  floating  frame  of  reference 
formulation.  The  use  of  such  a  local  element  coordinate  system  leads  to  a  highly  nonlinear 
expression  for  the  inertia  forces  as  the  result  of  the  large  element  rotation.  In  this  inves¬ 
tigation,  an  absolute  nodal  coordinate  formulation  is  presented  for  the  large  rotation  and 
deformation  analysis  of  three  dimensional  beam  elements.  This  formulation  leads  to  a  con¬ 
stant  mass  matrix,  and  as  a  result,  the  vectors  of  the  centrifugal  and  Coriolis  forces  are 
identically  equal  to  zero.  The  formulation  presented  in  this  paper  takes  into  account  the 
effect  of  rotary  inertia,  torsion  and  shear  effects,  and  ensures  continuity  of  the  slopes  as  well 
as  the  rotation  of  the  beam  cross  section  at  the  nodal  points.  Using  the  proposed  formulation 
curved  beams  can  be  systematically  modeled. 
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1  INTRODUCTION 


In  the  two  dimensional  analysis,  the  center  line  of  a  beam  element  can  be  used  to  completely 
describe  the  beam  kinematics  according  to  Euler-Bernoulli  beam  assumptions.  A  vector  that 
defines  the  location  of  an  arbitrary  point  on  the  beam  cross  section  in  terms  of  the  spatial 
longitudinal  coordinate  is  sufficient  for  the  determination  of  the  orientation  and  the  position 
vector  of  the  origin  of  a  Frenet  frame  [5,  12],  which  has  one  of  its  axes  tangent  to  the  center 
line,  and  the  other  two  axes  are  perpendicular  to  the  center  line.  The  effect  of  the  rotary 
inertia  of  the  beam  cross  section  can  be  systematically  accounted  for  using  Frenet  frame 
whose  orientation  is  completely  defined  by  the  derivatives  of  the  displacement  vector.  This 
was  the  basis  for  developing  the  absolute  nodal  coordinate  formulation  [11,  16,  18,  20],  which 
does  not  require  the  interpolation  of  finite  rotation  coordinates  as  it  is  the  case  in  some  finite 
element  procedures.  It  was  also  shown  that  the  absolute  nodal  coordinate  formulation  can 
account  for  the  rotary  inertia  effect  and  at  the  same  time  the  mass  matrix  remains  constant 
[16,  23].  The  fact  that  the  mass  matrix  remains  constant  becomes  crucial  in  developing  an 
efficient  algorithm  for  solving  the  multibody  dynamic  equations.  Using  this  property,  an 
optimum  sparse  matrix  structure  can  be  obtained  using  Cholesky  coordinates  [19,  24]. 

The  problem  of  three  dimensional  beams  requires  more  careful  consideration  since  the 
motion  of  the  beam,  even  within  the  Euler-Bernoulli  beam  assumptions,  can  not  be  com¬ 
pletely  described  by  the  displacement  of  its  center  line.  The  center  line  represents  a  spatial 
curve,  as  shown  in  Fig.  1.  The  location  of  an  arbitrary  point  on  this  spatial  curve  in  a 
coordinate  system  XYZ  is  defined  by  the  vector  r(s),  where  s  is  the  arc  length.  While 
the  vector  r(s)  can  be  used  to  define  a  Frenet  frame  which  has  three  orthogonal  vectors; 
tangent,  normal,  and  binormal,  such  a  representation  fails  to  capture  simple  beam  motion. 
For  instance,  the  rotation  of  the  beam  about  its  cross  section  as  a  rigid  body  can  not  be 
described  by  the  vector  r  since  this  vector  remains  constant  throughout  this  simple  motion. 
This  problem  has  been  avoided  in  the  finite  element  and  multibody  literature  by  introducing 
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a  local  coordinate  system.  In  the  finite  element  literature,  convected  coordinate  systems  are 
used  for  the  finite  element  [1,  15].  The  deformation  of  the  element  can  be  defined  in  the 
convected  system  whose  orientation  can  be  described  using  three  independent  parameters. 
In  the  multibody  literature,  the  floating  frame  of  reference  formulation  is  used  [2,  6,  9,  17]. 
In  this  formulation,  a  coordinate  system  is  introduced  for  the  flexible  body.  This  coordinate 
system  can  be  used  to  describe  the  gross  body  motion.  The  body  deformation  is  defined  in 
the  body  coordinate  system  which  can  capture  the  rotation  of  the  beam  element  about  its 
own  axis.  Another  benefit  of  using  the  floating  frame  of  reference  formulation  is  the  exact 
representation  of  the  rigid  body  motion  even  when  conventional  non-isoparametric  elements 
such  as  beams  and  plates  are  used.  With  this  formulation,  the  finite  element  has  zero  strain 
under  an  arbitrary  rigid  body  motion. 

It  is  important  to  point  out  that  most  of  existing  finite  element  formulations,  including 
large  rotation  vector  formulations  [3,  21,  22],  lead  to  a  highly  nonlinear  mass  matrix  when 
large  rotation  and  deformation  of  three  dimensional  beam  elements  are  considered.  It  is  the 
objective  of  this  investigation  to  develop  an  absolute  nodal  coordinate  formulation  for  three 
dimensional  beams  that  undergo  large  rotations  and  deformations.  It  is  shown  that  the  mass 
matrix  of  the  element  is  constant,  and  as  a  consequence,  the  centrifugal  and  Coriolis  inertia 
forces  are  identically  equal  to  zero.  The  effects  of  the  rotary  inertia,  torsion  and  shear  are 
automatically  accounted  for.  The  formulation  of  the  mass  matrix  and  the  elastic  forces  is 
presented  and  the  choice  of  the  element  nodal  coordinates  is  discussed. 


2  DISPLACEMENT  FIELD 

The  order  of  the  polynomial  used  in  the  assumed  displacement  fields  in  the  finite  element 
analysis  depends  on  the  expected  shape  of  the  element  deformation.  In  the  conventional 
finite  element  analysis,  a  linear  displacement  field  is  assumed  to  interpolate  the  longitudinal 
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deformation,  while  a  cubic  polynomial  is  used  for  the  transverse  deflections  of  the  beam 
element.  This  kind  of  assumed  displacement  field  is  suitable  when  using  a  beam  local  coor¬ 
dinate  system  as  in  the  case  of  the  incremental  procedure  or  the  finite  element  floating  frame 
of  reference  formulation  that  is  widely  used  in  flexible  multibody  simulation.  However,  in 
the  non-incremental  absolute  nodal  coordinate  formulation,  the  concept  is  different  since  the 
chosen  displacement  field  represents  both  rigid  and  flexible  body  motions.  This  displacement 
field  is  defined  in  the  global  system  and  accounts  for  the  coupling  between  the  rigid  body 
motion  and  the  elastic  deformation.  For  a  three  dimensional  beam  element,  we  assume  a 
displacement  field  using  the  following  polynomials: 
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n  - 

ao  +  aix  +  a2y  +  a^z  -1-  a^xy  +  a^xz  Y]  akX^~^ 

ri 

fc=6 

jl 

r2 

bo  -f-  bix  -f  b2y  +  b^z  +  b^xy  H-  b^xz  -1-  ^  bkx'^~‘^ 

/c=6 

n 

rs 

Co  -f  cix  +  C2y  +  csz  +  c^xy  +  c^xz  +  ^  Ckx'^~‘^ 

k=e 

where  r  is  the  global  position  vector  of  an  arbitrary  point  on  the  beam  cross  section,  Uj,  6j, 
and  Ci  are  the  polynomial  coefficients,  and  x,  y  and  z  are  the  spatial  coordinates  defined  in 
a  chosen  beam  coordinate  system.  Here  x  is  assumed  to  be  the  spatial  coordinate  along 
the  beam  axis  {0  ^  x  ^  1),  where  I  is  the  length  of  the  beam  element.  The  order  of  the 
polynomials  in  the  preceding  equation  can  be  chosen  depending  on  the  magnitude  of  the 
deformation  expected  from  the  element.  Note  that  the  polynomials  used  to  interpolate  the 
three  components  of  the  displacements  have  the  same  order  since  the  vector  r  is  defined  in 
the  global  coordinate  system.  Furthermore,  in  order  to  account  for  the  rotary  inertia  and 
shear  effects,  the  displacement  field  is  assumed  to  depend  on  y  and  z.  Since  the  cross  section 
dimensions  of  the  beam  element  are  assumed  to  be  small  compared  to  the  element  length, 
the  displacement  field  is  assumed  to  depend  only  linearly  on  the  spatial  coordinates  y  and 
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z.  The  preceding  equation  can  also  be  written  in  a  matrix  form  as  follows;  '  ■ 

(2) 

where  a,  b,  and  c  are  the  vectors  of  the  polynomial  coefficients,  and  Si  is  a  row  vector  that 
defines  the  space-dependent  coefficients  of  the  polynomials  of  Eq.  (1). 

The  order  of  the  interpolating  polynomials  used  for  the  assumed  displacement  field  de¬ 
pends  on  the  chosen  number  of  the  beam  nodal  coordinates.  The  choice  of  the  nodal  coor¬ 
dinates  is  as  important  as  the  choice  of  the  displacement  fields  since  it  contributes  to  the 
definition  of  the  element  stiffness  and  inertia  forces.  A  combination  of  global  position  and 
slope  coordinates  can  always  be  chosen  as  nodal  coordinates  as  will  be  demonstrated  in 
details  in  later  sections  of  this  paper. 
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3  DEFINITION  OF  BEAM  CROSS  SECTION 


In  Euler-Bernoulli  beam  theory,  the  cross  section  is  assumed  to  remain  rigid  and  perpendic¬ 


ular  to  the  center  line.  Therefore,  the  normal  to  the  cross  section  is  defined  by  the  tangent 
vector  dv/dx.  In  the  assumed  displacement  field  defined  in  the  preceding  section,  the  beam 
cross  section  does  not  remain  perpendicular  to  the  center  line,  and  therefore,  the  cross  sec¬ 


tion  is  not  defined  by  the  vector  dr/dx  tangent  to  the  center  line.  Using  this  displacement 

dv  dv 

field,  it  can  be  shown  that  for  a  given  x,  —  and  —  are  independent  of  y  and  z.  It  can  also 

oy  oz 

dv  dv 

be  shown  that  —  and  —  are  two  independent  vectors  (not  necessarily  orthogonal)  that 
oy  oz 

define  the  cross  section  of  the  beam  element.  To  this  end,  an  arbitrary  vector  r^.  is  defined  in 


the  cross  section  as  shown  in  Fig.  2.  Using  the  displacement  field  defined  in  the  preceding 
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section,  it  can  be  shown  that 


(I2  “h  XCI4 

as  +  xa^ 

=  rp  -  rp  = 

62  +  X64 

y  + 

63  +  xb^ 

C2  +  XC4, 

Cs  +  XC5 

where  rp  is  the  global  position  vector  of  an  arbitrary  point  P  on  the  cross  section  with 
coordinates  {x,y,z),  and  rp  is  the  global  position  vector  of  a  corresponding  point  P  with 
coordinates  {x,  0, 0)  on  the  center  line  of  the  beam. 


It  can  be  shown  that 


a2  +  xci4^ 

03  +  xas 

dr 

dr 

dy 

62  +  xb/^ 

’  dz  ” 

bz  +  xbs 

C2  +  XC4 

C3  +  XC5 

The  preceding  two  equations  (Eqs.  3  and  4)  show  that 


dr  dr 


(4) 


which  shows  that  an  arbitrary  vector  drawn  on  the  cross  section  can  be  expressed  as  a  linear 


combination  of  the  two  vectors  7^  and 

oy 


dr 

dz 


The  fact  that 


^  and  ^  are  independent 
dy  dz 


of  y  and  ^  can  be  used,  with  the  help  of  Gram-Shmidt  orthogonalization  process  [12],  to 


determine  the  following  two  orthonormal  vectors  on  the  cross  section: 


n*  = 


b.= 


Fz  -  hry 
Fz 


(5) 


where 


(FzfF^ 

ri/^F. 


(6) 


dv 

n^  is  a  unit  vector  along  r.y  =  — ,  while  bj  is  a  unit  vector  which  is  a  linear  combination  of 

dy 

dv 

the  vectors  =  —  and  n^.  In  order  to  complete  the  orthogonal  triad  for  the  cross  section, 
dz 

a  normal  to  the  beam  cross  section  can  be  obtained  by  using  the  cross  product  of  the  vectors 


n.5  and  bs  as  follows: 


tg  X  bj 


(7) 
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The  vectors  and  bj  constitute  a  right-handed  orthonormal  set,  which  completely 

defines  the  orientation  of  the  beam  cross  section. 


4  KINEMATICS  OF  THE  CROSS  SECTION 

As  described  in  the  preceding  section,  the  three  orthogonal  vectors  and  can  be 

systematically  defined  in  terms  of  the  vectors  Vy  and  obtained  using  the  displacement 
field  introduced  in  Section  2.  In  the  case  of  a  rigid  body  motion  of  the  finite  element  it  can 
be  shown  that  and  are  two  orthogonal  unit  vectors,  and  as  a  consequence,  h  in  Eq.  (6) 
is  identically  equal  to  zero,  and 


lis  =  Ty,  and  hg  =  r^.  (8) 

Furthermore,  in  the  case  of  rigid  body  motion  of  a  straight  beam,  remains  in  the  direction 
of  the  center  line  of  the  beam.  Therefore,  in  the  case  of  a  general  displacement  that  includes 
deformations,  the  difference  between  n*  and  r^,  and  bj  and  is  mainly  due  to  the  beam 
deformation.  In  fact  h  in  Eq.  (6)  can  be  used  as  a  measure  of  this  difference  since  h  is 
identically  equal  to  zero  under  an  arbitrary  rigid  body  motion.  If  the  deformation  within 
the  element  is  small,  then  it  is  reasonable  to  assume  that 

/i  ~  0, 


and  in  this  case 


n^.  Si  Tj,,  and  b^,  ss  r^. 


Using  this  assumption,  it  follows  from  the  definition  of  presented  in  the  preceding  section 
that 

Ir^l  «  constant, 


which  implies  that  the  cross  section  remains  rigid. 
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It  is  important  to  point  out  that  the  assumption  used  in  this  section  is  also  commonly 
used  in  the  finite  element  literature  for  beam  elements  that  account  for  the  effect  of  rotary 
inertia  and  torsion  [14]  as  will  be  discussed  in  Section  9. 


MEASURES  OF  TORSION  AND  SHEARS 


In  Euler-Bernoulli  beam  theory,  it  is  assumed  that  the  cross  section  of  the  beam  remains 
perpendicular  to  the  beam  center  line,  and  as  a  result,  the  shear  deformation  effect  is  ne¬ 
glected.  In  this  case,  the  vector  normal  to  the  cross  section  coincides  with  the  tangent  to 
the  center  line  defined  by  Serret-Frenet  equations. 

Serret-Frenet  frame  is  a  coordinate  system  which  defines  three  vectors;  tangent,  normal, 
and  binormal  to  the  beam  center  line.  The  normal  and  binormal  vectors  define  the  so  called 
normal  plane  [12]  which  is  normal  to  the  center  line.  The  normal  plane  of  Serret-Frenet 
frame  is,  therefore,  the  beam  cross  section  if  Euler-Bernoulli  beam  theory  is  used.  The 
Serret-Frenet  frame  is  defined  by  three  orthogonal  unit  vectors;  t,  n,  and  b.  The  unit  vector 
t  is  tangent  to  the  beam  center  line  and  is  defined  as 

dr 


t  = 


ds' 


(9) 


where  ds  is  the  arc  length  of  an  infinitesimal  segment  of  the  center  line  defined  as 


ds  =  \l  {r^f  r^dx, 


(10) 


dr 


where  Thus,  Eq.  (9)  can  be  written  as 

kJ  Jb 


t 


Yxf 

The  unit  vector  n  is  normal  to  the  beam  center  line  and  is  defined  as 

I  dt  1  dt  dx 


(11) 


n  = 


Kds  ndx  ds’ 


(12) 
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The  rate  at  which  the  curve  twists  out  of  its  osculating  plane  is  called  the  torsion  r,  which 
is  defined  as  [12] 

cfb 
ds 

and  it  represents  the  twisting  shear  deformation  of  the  beam.  However,  as  pointed  out  in  the 
introduction  of  this  paper,  the  Serret-Frenet  description  fails  to  capture  the  rotation  of  the 
cross  section  of  a  straight  beam  about  its  own  axis  since  such  a  rotation  does  not  contribute 
to  the  change  of  the  vector  r. 

As  already  mentioned,  the  plane  defined  by  {n,  b}  is  perpendicular  to  the  center  line 
of  the  beam,  while  the  plane  defined  by  {ns,b^}  represents  the  cross  section  of  the  beam 
defined  by  the  assumed  displacement  field  presented  in  Section  2.  In  the  undeformed  beam 
configuration  (rigid  body  motion),  or  in  the  case  of  no  shear  effect,  the  two  planes  are  parallel. 
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In  the  case  of  a  general  displacement,  the  two  planes  are  not  parallel  due  shear  deformations. 
The  scalar 

V  =  t-  ts=  t-  (n^  X  bs)  =  det  nj  ,  (18) 

gives  an  idea  about  the  shear  deformation  of  the  beam  as  shown  in  Fig.  3.  For  instance,  if 
V  =  1,  there  is  no  shear  deformation,  and  if  it  is  not  equal  to  one,  the  amount  of  shear  is 
inversely  proportional  tor;(— l^n^l).  Figure  4  illustrates  a  general  deformation  of  the 
beam  by  assuming  arbitrary  values  for  the  polynomial  coefficients.  Also  shown,  in  Fig.  4,  is 
the  angle  7  which  is  the  arc  cosine  of  v  at  two  different  points  on  the  beam  =  x/l).  Later 
in  this  paper,  the  strain  energy  due  to  the  shear  effect  will  be  discussed. 

Using  the  assumed  displacement  field  of  the  beam  presented  in  Section  2,  the  effect  of 
the  torsion  can  be  easily  accounted  for  using  the  two  unit  vectors  Uj  and  bs  of  the  beam 
cross  section.  As  a  measure  of  torsion  for  a  short  beam,  one  can  use  the  rotation  of  the  cross 
section  at  an  arbitrary  point  on  the  beam  center  line  with  respect  to  the  reference  plane  at 
node  A.  There  exists  an  angle  Px  defined  by  the  dot  product 

cos Px  =  ris  •  i^sA,  smpx  =  ns-hsA,  (19) 

where  and  bj^  are  the  two  unit  vectors  that  define  the  reference  plane.  Note  that  Px 
is  equal  to  zero  in  the  case  of  rigid  body  motion,  and  as  a  consequence,  it  can  be  used  as  a 
measure  for  the  torsion. 

6  BEAM  INERTIA 

In  the  three  dimensional  analysis  of  the  large  rotation  problem,  both  the  finite  element 

incremental  approach  and  large  rotation  vector  formulation  lead  to  a  highly  nonlinear  mass 
matrix.  The  absolute  nodal  coordinate  formulation  on  the  other  hand  leads  to  a  constant 
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mass  matrix  and  automatically  accounts  for  the  shear  and  torsional  effects. 

The  mass  matrix  of  the  three  dimensional  beam  element  can  be  obtained  using  the 
absolute  nodal  coordinate  formulation  and  the  following  expression  of  the  kinetic  energy': 

T  =  i  /  (20) 

V 

where  r  is  the  global  position  vector  of  an  arbitrary  point,  and  p  and  V  are  respectively  the 
mass  density  and  volume  of  the  beam  element.  The  vector  r  is  defined  in  Section  2  by  Eq. 
(1)  using  the  coefficients  of  the  interpolation  polynomials.  These  coefficients  can  be  replaced 
by  global  coordinates  and  slopes  at  the  nodes,  as  explained  in  later  sections  of  this  paper. 
In  this  case,  the  global  position  vector  of  the  arbitrary  point  can  be  expressed  in  terms  of  a 
vector  of  element  nodal  coordinates  e  and  the  shape  function  S  as 

r  =  Se,  (21) 


It  follows  that 


f  =  Se. 


(22) 


By  substituting  Eq.  (22)  into  Eq.  (20),  one  obtains 


e, 


(23) 


which  is  a  simple  quadratic  form  in  the  velocities.  Thus  the  element  mass  matrix  is  defined 
as 


M  = 


(24) 


V 

The  above  integration  defines  a  constant  mass  matrix  which  only  depends  on  the  inertia 
properties  and  dimensions  of  the  beam.  Using  the  fact  that  the  mass  matrix  obtained  using 
the  absolute  nodal  coordinate  formulation  is  constant,  efficient  numerical  procedures  can  be 
used  to  obtain  an  optimum  sparse  matrix  structure  for  the  resulting  multibody  dynamic 
equations  [19,  24] 


10 


The  inclusion  of  the  effect  of  the  rotary  inertia  of  the  three  dimensional  formulation 
presented  in  this  paper  can  be  demonstrated  by  using  Eq.  (3).  To  this  end,  we  write 

Tp=  Tp  + is,  (25) 

where  P  is  an  arbitrary  point  on  the  beam  center  line.  Using  the  preceding  equation,  the 
beam  kinetic  energy  can  be  written  as 

T  =  P{if  +  iJ){tp  +  i,)dV  (26) 

p  [2rprs  +  ij f  j]  dxdydz 

Ki  V 

The  first  term  in  the  above  integration  is  the  mass  matrix  in  the  case  in  which  the  beam 
rotary  inertia  is  neglected,  while  the  second  term  accounts  for  the  effect  of  the  rotary  inertia. 


7  ELASTIC  FORCES 

If  the  deformation  within  the  beam  is  large,  the  expression  for  the  nonlinear  strain-displacement 
relationship  must  be  used  to  formulate  the  elastic  forces  of  the  beam  element  [4]  in  the  non- 
incremental  absolute  nodal  coordinate  formulation.  However,  if  the  size  of  the  element  is 
chosen  to  be  small  such  that  the  deformation  within  the  element  remains  small,  the  linear 
strain-displacement  relationship  can  be  used  in  the  large  deformation  analysis  of  flexible  bod¬ 
ies  using  the  absolute  nodal  coordinate  formulation  [20,  24].  In  this  section,  as  an  example 
for  formulating  the  elastic  forces,  the  case  of  small  element  deformation  is  considered. 

In  the  finite  element  analysis,  the  strain  energy  of  a  beam  is  defined,  in  general,  in  terms 
of  six  components;  one  axial  force,  two  bending  moments,  two  shear  forces,  and  one  torsional 
moment  [14].  However,  in  Euler-Bernoulli  beam  theory,  the  shear  forces  are  not  considered 
assuming  that  the  cross  section  plane  remains  perpendicular  to  the  beam  centerline  after 
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deformation.  Thus,  the  strain  energy  in  the  case  of  small  deformation  can  be  defined  as  [10] 

EA  + Ei„  {py + Eh.  [^y+ 

Gk0^  +  Ok0^  +  G/ix  {^7) 
where  Ux,  Uy,  and  are  respectively  the  x—,  y—,  and  z— component  of  the  beam  deflection, 
Px,  Py,  and  Pz  are  the  shear  angles,  k  is  the  Timoshenko  shear  factor  [8,  10],  E  and  G  are 
respectively  the  moduli  of  elasticity  and  rigidity,  Ixx,  lyy,  and  Ez  are  the  second  moments  of 
area,  and  A  is  the  beam  cross  section  area.  In  order  to  evaluate  the  beam  strain  energy,  a 
local  beam  frame  can  be  used.  However,  it  is  important  to  point  out  that  the  beam  strain 
energy  can  be  also  evaluated  using  the  inertial  coordinate  system  directly  instead  of  using  a 
local  frame  by  utilizing  continuum  mechanics  theories  [4,  21,  22,  23]. 

Previously,  two  coordinate  systems  were  introduced  in  this  paper;  the  cross  section  frame 
and  the  Serret-Frenet  frame.  In  order  to  define  the  beam  longitudinal  and  transverse  deflec¬ 
tions,  a  vector  d,  shown  in  Fig.  5,  is  defined  as 

d  =  r/s  -  r^,  (28) 

where  rp  is  the  position  vector  of  an  arbitrary  point  P  on  the  beam  centerline,  and  is 
the  position  vector  of  node  A  as  shown  in  Fig.  5.  Considering  Eq.  (21),  Eq.  (28)  can  be 
re-written  as 

d-(S(a:,0,0)-S(0,0,0))e.  (29) 

As  shown  in  Fig.  5,  the  vector  d  represents  the  location  of  a  point  P  on  the  beam  center 
line  with  respect  to  node  A,  the  deflection  components  can  be  defined  as  follows: 


'^x 

X 

u  — 

Uy 

=  A^d- 

0 

Uz 

0 

where  x  is  the  axial  coordinate  of  point  P  in  the  undeformed  state,  A  is  a  transformation 
matrix  which  can  be  expressed  in  terms  of  the  unit  vectors  t^,  n*,  and  b*  defined  at  node  A. 


I  ( 
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It  follows  from  the  preceding  equation  that 


Ux  ~  ‘  d  X 

Uy  -  ’  d  ^  • 

Uz  —  bsA  ■  d 


Equation  (31)  defines  the  axial  and  transverse  components  of  the  beam  defiection  which  are 
required  for  evaluating  the  first  three  terms  of  the  strain  energy. 

There  are  different  approaches  that  can  be  used  to  measure  the  torsion  and  shear.  Figures 
6a,  b,  and  c  show  a  schematic  diagram  for  each  of  the  three  rotation  angles  of  the  beam  cross 
section.  Using  the  difference  in  orientation  between  the  cross  section  frame  and  Serret-Frenet 
frame  depicted  in  these  figures,  one  can  conclude  the  following: 


cosPx  = 

sinpx  =  t^s-^sA 

where  px  is  the  rotation  of  the  cross  section  relative  to  a  reference  plane  which  is  assumed 
to  be  at  node  A,  and  are  the  unit  vectors  that  define  the  reference  plane  at  node  A. 
It  can  be  shown  using  Eq.  (32)  that 

(^)  "" 


where  n,^  = 


dus 

dx 


The  shear  angles  can  be  defined  as  the  rotations  of  the  beam  cross  section 


about  the  normal  and  binormal  vectors  of  Serret-Frenet  frame.  Considering  Fig.  6b,  the 


shear  angle  jSy  can  be  defined  as 


Py  =  sin  ^  {hs  ■  t) . 


(34) 


Similarly,  the  shear  angle  Pz,  according  to  Fig.  6c,  is  defined  as 

Pz  =  sin“^  (-ns  •  t) . 


(35) 
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Using  Eq.  (27),  the  definition  of  the  longitudinal  and  transverse  deformations,  the  torsion, 
and  the  shear  angles;  the  strain  energy  of  the  beam  can  be  written  as 


where  e  is  the  vector  of  nodal  coordinates,  K  is  the  stiffness  matrix.  In  the  current  analysis, 
the  beam  stiffness  matrix  is  highly  nonlinear  in  the  nodal  coordinates.  Differentiating  the 
strain  energy  once  with  respect  to  the  nodal  coordinates  will  lead  to  the  beam  elastic  forces 
Q  defined  as 


Q  — 


m 

de 


(37) 


In  order  to  ensure  that  the  proposed  element  meets  the  convergence  requirements,  an  eigen¬ 
value  test  has  to  be  performed  on  the  stiffness  matrix  [7].  The  stiffness  matrix  must  have 
the  exact  number  of  rigid  body  modes  (3  for  2-D  and  6  for  3-D  elements)  despite  the  fact 
that  the  stiffness  matrix  is  function  of  the  nodal  coordinates.  In  other  words,  the  stiffness 
matrix  should  include  exactly  6  zero  eigenvalues.  If  the  stiffness  matrix  has  zero  eigenvalues 
more  than  the  number  of  the  rigid  body  modes,  this  implies  that  the  element  has  zero- 
energy  deformation  modes  which  do  not  have  corresponding  restoring  forces  represented  in 
the  beam  strain  energy.  Such  an  element  will  not  converge  to  a  correct  solution.  Also  the 
stiffness  matrix  should  produce  zero  strain  energy  (zero  elastic  forces)  in  the  case  of  rigid 
body  motion. 


8  NODAL  COORDINATES 

In  the  non-incremental  absolute  nodal  coordinate  formulation,  the  displacement  field  of  the 
finite  element  is  expressed  in  terms  of  a  set  of  nodal  coordinates  that  consists  of  global 
displacement  and  slope  coordinates.  The  number  of  these  coordinates  depends  on  the  order 
of  the  polynomials  used.  The  displacement  field  presented  in  Section  2  is  assumed  to  be 
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linear  in  y  and  2:  and  it  can  assume  any  order  in  x.  For  example,  if  cubic  polynomials  in 
X  are  used,  one  needs  24  nodal  coordinates  for  the  beam  element.  This  number  can  be 
reduced  by  reducing  the  order  of  the  polynomials.  However,  by  increasing  the  order  of  the 
polynomials,  more  deformation  modes  within  one  element  can  be  obtained.  The  use  of  such 
a  higher  order  element  may  lead  to  a  significant  reduction  in  the  number  of  finite  elements 
required  to  model  a  structure.  Figure  7  shows  deformed  shapes  which  can  be  obtained  using 
one  element  with  polynomials  cubic  in  x.  The  shape  in  Fig.  7a  is  obtained  using  the  following 
polynomial  coefficients: 

af  =  1  1  0  0  0  0  -0.27  0.018  ,  bf  =  1  1  0  0  0  0  -0.04  -0.004  , 

cf  =  0  1  0  0  0  0  -0.11  0.004  , 

while  the  shape  presented  in  Fig.  7b  is  obtained  using  the  following  coefficients; 

a J  =  1  -3  0  0  0  0.77  0.048  ,  =  1  -2  0  0  0  0  -0.02  0.018  , 

C2  =  0  1  0  0  0  0  -0.13  0.002  . 

Such  deformed  shapes  can  not  be  obtained  with  one  or  two  elements  using  linear  interpola¬ 
tion.  Depending  on  the  order  of  the  polynomial  used  and  the  number  of  nodes  selected  for 
the  element,  one  can  choose  coordinates  using  the  following  global  variables: 

dr  dr  dr 

’  dx'  dy'  dz 

These  variables  represent  12  coordinates  at  a  given  point  on  the  center  fine  of  the  element. 
One  does  not  need  to  use  all  these  coordinates  at  a  given  node.  One  may  also  choose  to 
use  another  set  of  coordinates  which  include  a  linear  combination  of  the  global  slopes.  For 
example,  the  following  three  variables  can  be  used 

_  dr2  dri 

^  dy  dz'  ^  dz  dx'  ^  dx  dy 
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One  can  show  that  in  the  case  of  rigid  body  motion  'dx,'dy,  and  can  be  related  to  the 
known  orientation  parameters  used  in  the  rigid  body  dynamics.  For  example,  'dx,'dy,  and  dz 
can  be  expressed  in  terms  of  Euler  parameters,  6q,  61,82,  and  9z,  in  the  case  of  rigid  body 
motion  as  follows: 


=4^0^1,  =  400^2,  ^9z=4M3, 

and  they  are  related  to  Rodriguez  parameters  as  follows  [18]: 

471  .0  472  ,  473 


1  +  7"’ 


1+7^’ 


•dz 


1  +  7^’ 


where  71,72,  and  73  are  Rodriguez  parameters  and  7  =  a/ti  +  tI  +  73  • 

A  large  rotation  vector  can  be  defined  as  v  sin  6  where  v  is  a  unit  vector  along  the  axis 
of  rotation  and  6  is  the  angle  of  rotation.  It  can  be  shown  that,  in  the  case  of  rigid  body 

lT 


motion,  the  set 'd 


dx  'dy  tDz 


can  be  expressed  in  terms  of  the  rotation  vector  as 


■d  =2vsin0. 


9  BASIC  ASSUMPTIONS 

In  several  finite  element  formulations  for  the  large  rotation  and  deformation  analysis  of 
beams,  the  equations  of  motion  are  developed  by  assuming  that  the  beam  cross  section 
remains  rigid  [21,  22],  while  in  reality,  the  beam  cross  section  does  not  remain  rigid  when 
the  beam  deforms.  The  use  of  the  assumption  of  the  rigidity  of  the  cross  section  introduces 
difficulties  when  these  methods  are  generalized  to  the  case  of  plates  and  shells  where  the 
rigidity  assumption  can  not  be  used  and  is  no  longer  valid.  The  non- incremental  absolute 
nodal  coordinate  formulation  relaxes  this  assumption,  and  therefore,  it  can  be  easily  extended 
to  study  plate  and  shell  problems. 

It  was  shown  in  Section  4  that  the  length  of  a  vector  on  the  cross  section  of  the  beam 
element  remains  constant  under  an  arbitrary  rigid  body  displacement.  In  the  general  case 
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of  arbitrary  displacements  and  deformation,  the  change  in  the  length  of  this  vector  is  small 
and  is  due  to  the  deformation.  It  is  the  objective  of  this  section  to  demonstrate  that  this 
hypothesis  is  commonly  used  in  the  finite  element  literature.  To  this  end,  we  use  the  con¬ 
ventional  finite  element  shape  function  used  for  beams  and  for  simplicity  we  consider  the 
case  of  two-dimensional  beam  element.  Nonetheless,  the  main  conclusions  obtained  using 
the  simpler  two  dimensional  beam  model  apply  to  the  conventional  three  dimensional  beam 
element  used  in  the  incremental  procedure. 

In  the  assumed  displacement  field  of  the  conventional  beam  element,  six  nodal  coordinates 
are  used  for  the  two-node  beam  element  (3  coordinates/node;  two  displacement  coordinates, 
and  one  infinitesimal  rotation  coordinate).  The  element  shape  function  of  the  Euler-Bernoulli 
beam  element  is  given  by 

^0  0 
0  3(0" -2  (O'  «[(0"-(0'] 

(38) 


s  = 


^-10  0 
0  1-3(0" +  2  (O'  2(0'  + (O'] 


where  ^  =  x//,  and  the  vector  of  nodal  coordinates  is 


Prs  Po  P  P, 


where 


ei  = 

^2  —  ’’21^=0  > 

63  = 

dr-i 

dx 

64  = 

^5  ^2 1^=1  ! 

66 

dr2 

dx 

Therefore,  the  coordinates  of  an  arbitrary  point  on  the  center  line  of  the  beam  element 
defined  in  a  convected  coordinate  system  is 


Note  that  this  element  does  not  account  for  the  effect  of  the  rotary  inertia.  In  the  literature, 
the  effect  of  the  rotary  inertia  is  considered  by  modifying  the  assumed  displacement  field  as 
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follows: 


r  =  Se+ 


y- 


(39) 


—dr2ldx 
drxfdx 

Since  the  shape  function  of  Eq.  (38)  is  independent  of  y,  the  vector  tangent  to  the  center  line 
is  defined  as  dv/dx  which  can  be  used  to  define  the  normal  vector  n  =  —dr2ldx  dvxjdx 
that  appears  in  the  preceding  equation.  If  the  element  length  and  the  axial  deformations  are 
assumed  to  be  small,  then 

1.  (40) 


dx 

It  follows  that  a  vector  on  the  cross  section  is  defined  in  the  convected  coordinates  system 
as 

—dr2/dx 
dvi/dx 

which  implies  that  |rs|  =constant  if  and  only  if  n  is  a  unit  vector.  This  is  not,  in  general, 
the  case  since  Eq.  (40)  is  not  strictly  imposed.  To  strictly  impose  the  condition  of  Eq.  (40), 
one  must  write 


n  = 


—dr2/dx 

dri/dx 


{dri/dx)^  +  {dr2fdxf 

which  will  produce  a  shape  matrix  which  is  nonlinear  in  the  element  nodal  coordinates. 


10  SUMMARY  AND  CONCLUSIONS 

A  non-incremental  absolute  nodal  coordinate  formulation  for  three  dimensional  beam  ele¬ 
ments  is  presented  in  this  paper.  The  element  displacement  field  is  assumed  to  be  linear  in 
y  and  z  coordinates  of  the  cross  section  of  the  beam,  while  it  can  assume  higher  order  in  x. 
It  can  be  shown  that  this  assumed  displacement  field  can  describe  exact  rigid  body  motion, 
and  as  a  consequence,  it  leads  to  zero  strain  energy  under  an  arbitrary  rigid  body  displace¬ 
ment.  The  kinematics  of  the  element  cross  section  is  thoroughly  examined  and  compared 
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with  the  Serret-Frenet  frame.  It  is  shown  that  the  tangent  of  the  center  line  of  the  element 
does  not  remain  perpendicular  to  the  element  cross  section,  thereby  demonstrating  that  the 
new  element  can  account  for  the  shear  deformation  effect  as  well  as  the  torsion.  The  shear 
and  torsion  angles  are  defined  in  terms  of  the  deviation  of  the  element  cross  section  frame 
from  the  Serret-Frenet  frame.  The  formulation  of  the  inertia  and  elastic  forces  of  the  finite 
element  using  the  absolute  nodal  coordinate  formulation  was  discussed  in  the  paper.  It  is 
shown  that  it  is  feasible  to  obtain  an  element  that  has  a  constant  mass  matrix  and  at  the 
same  time  accounts  for  the  rotary  inertia,  shear,  and  torsion  effects.  This  property  is  an 
important  feature  of  the  absolute  nodal  coordinate  formulation  since  most  existing  methods 
used  for  the  nonlinear  large  rotation  analysis,  including  the  incremental  methods,  large  ro¬ 
tation  vector  formulations,  and  the  floating  frame  of  reference  formulation,  lead  to  a  highly 
nonlinear  mass  matrix  when  three-dimensional  beams  are  considered. 

Global  slopes  or  linear  combination  of  these  slopes,  in  addition  to  displacement  coor¬ 
dinates,  can  be  used  as  the  element  nodal  coordinates.  The  number  of  the  element  nodal 
coordinates  depends  on  the  order  of  the  interpolating  polynomials  presented  in  Section  2. 
The  assumptions  used  in  developing  the  three  dimensional  beam  element  presented  in  this 
paper  are  also  discussed  and  it  is  demonstrated  that  these  assumptions  are  consistent  with 
the  assumptions  that  have  been  used  in  other  methods  that  employ  the  conventional  three 
dimensional  beam  elements.  Some  other  important  features  of  the  method  proposed  in  this 
paper  can  be  summarized  as  follows: 

1.  The  method,  in  general,  relaxes  the  assumption  of  the  rigidity  of  the  cross  section  of 
the  beam  element,  and  therefore,  it  can  be  extended  to  the  analysis  of  plates  and  shells. 

2.  The  method  leads  to  isoparametric  elements  that  can  be  easily  used  in  the  analysis  of 
curved  structures. 

3.  The  method  does  not  require  the  interpolation  of  finite  rotation  coordinates.  It  is 
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known  that 


(a)  for  a  given  configuration  and  a  given  set  of  orientation  coordinates,  different 
sequences  with  different  values  for  the  orientation  coordinates  can  be  used  to 
describe  the  same  configuration. 

(b)  the  relationship  between  different  sets  of  orientation  coordinates  is  highly  nonlin¬ 
ear,  and  therefore,  a  linear  interpolation  of  one  set  of  finite  rotation  coordinates 
does  not  imply  the  same  order  of  approximation  for  another  set.  This  raises  some 
questions  with  regard  to  the  physics  used  in  the  interpolation  of  the  finite  rotation 
coordinates.  On  the  other  hand,  in  the  absolute  nodal  coordinate  formulation, 
only  the  shape  of  the  element  is  interpolated. 

4.  The  proposed  method  ensures  the  continuity  of  the  slopes  and  the  rotation  of  the  cross 
section.  As  a  consequence,  the  configuration  shown  in  Fig.  8  which  can  result  from  the 
interpolation  of  rotation  coordinates  only,  is  avoided.  This  configuration  is  the  result 
of  imposing  continuity  conditions  on  the  finite  rotations,  while  no  such  conditions  are 
imposed  on  the  global  slopes. 
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Fig.  3.  Coordinate  systems. 
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Fig.  4.  Rotation  of  the  beam  cross  section. 


(a)  Rotation  of  beam  cross  section  about  its  normal. 


(b)  Rotation  of  beam  cross  section  due  to  shear 


(c)  Rotation  of  beam  cross  section  due  to  shear 


Fig.  6.  Definition  of  torsion  and  shear. 
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Fig.  8.  Discontinuity  of  the  slopes. 


